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Abstract 

We study the asymptotic scaling properties of standard domain wall networks in several cosmo- 
logical epochs. We carry out the largest field theory simulations achieved to date, with simula- 
tion boxes of size 2048 3 , and confirm that a scale-invariant evolution of the network is indeed 
the attractor solution. The simulations are also used to obtain an accurate calibration for the 
velocity-dependent one-scale model for domain walls: we numerically determine the two free 
model parameters to have the values c K = 0.34 + 0.16 and k w = 0.98 + 0.07, which are higher 
precision than (but in agreement with) earlier estimates. 
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1. Introduction 

A key consequence of cosmological phase transitions is the formation of topological defects 
(HQ]. While cosmic strings have attracted most of the community's attention, domain walls are 
useful as a testbed case with which one can gather information relevant for other more complex 
defects (despite being tightly constrained by observations JHQl). Here we take advantage of ever- 
improving computing resources to carry out a large set of 204 8 3 high-resolution simulations of 
domain walls, using the standard Press-Ryden-Spergel (PRS) algorithm ||5|]. This is a follow-up 
on |@], where the results of simulations of size up to 1024 3 were presented, and we confirm and 
expand their results. 

Early generations of domain wall simulations 

earn a em in in 

found some hints 

for late-time deviations from the scale-invariant evolution, which would be the expected behav- 



ior 1114 11211 . Our previous work [6] found no such deviations, which provided support for the 
hypothesis that the earlier results were simply a consequence of the limited dynamical range of 
numerical simulations. We believe that the present work clearly confirms this. 

Macroscopic properties of defect networks can be accurately described by an analytic velocity- 
dependent model, first derived for cosmic strings lfl5l [l6 17]. The large-scale features of the 
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network are described by a characteristic scale L (which one can interchangeably think of as 
a typical defect separation or correlation length) and a microscopically averaged (root-mean- 
squared) velocity v. This has the advantages of tractability and conceptual simplicity but must 
include phenomenological parameters which parametrize our ignorance about certain dynamical 
mechanisms. The only way to accurately determine the correct values of these parameters is 
by employing large-scale numerical simulations to calibrate them. The main goal of the current 
work is precisely to improve this calibration. 



2. Numerical simulations 

We will study simple (single-field) domain wall networks in fiat homogeneous and isotropic 
Friedmann-Robertson-Walker (FRW) universes. (Throughout the paper we shall use fundamen- 
tal units, in which c — h — 1.) A scalar field (p with Lagrangian density 



\Y 



(1) 



provides the simplest case. By standard variational methods we obtain the field equation of 
motion (written in terms of physical time f) 

where V is the Laplacian in physical coordinates, H = a~ l (da/dt) is the Hubble parameter and a 
is the scale factor, which we generically assume to vary as a oc t 1 . In what follows we will study 
the network's evolution in several such cosmological epochs. 

We follow the procedure of Press, Ryden and Spergel [5], modifying the equations of motion 
in such a way that the thickness of the domain walls is fixed in co-moving coordinates. The 
reliability of this method has been numerically tested in previous work Q El EM] . In the PRS 
method, equation © becomes: 

d 2 <p ( d In a\d<f> 2 a^V 



+ a 



drj 1 \d\nrjj drj 



V 2 = -</—. (3) 



where 77 is the conformal time and a and /3 are constants: /3 — is used in order to have constant 
co-moving thickness and a = 3 is chosen to require that the momentum conservation law of the 
wall evolution in an expanding universe is maintained [5]. The specific parameters used in the 
simulations are 0o = 1, Vq = 7t 2 /2Wq, where Wq = 10 is the wall thickness; these choices are 
also justified by previous work on this algorithm iBlfll 181. 



Despite these previous results one may wonder whether the chosen wall thickness is suf- 
ficient to accurately calibrate the model. Figure [1] compares two series of ten 1024 3 mat- 
ter era runs with Wo = 10 and Wo = 100. Plotted are the ratios of the network velocities 
(R yv = (yv)M)=ioo/(y.)wo=io) and densities (R p = (p)wo=ioo/(p)wo=io)- After transients due to 
the choices of initial conditions, both of these ratios converge to unity, with statistical uncertain- 
ties well below 10%. This convergence is expected to be stronger wilth larger ensembles and/or 
larger boxes. 

Equation (0 is integrated using a standard finite-difference scheme. We assume the initial 
value of (j> to be a random variable between -<f>o and +(f>o and the initial value of dcp/dr] to be 

2 




Figure 1: The evolution of the velocity and density ratios (R yv and R p , defined in the text) for the average of two sets of 
ten 1024 3 matter era simulations with Wo = 100 and Wq = 10. As expected, both ratios become unity after a transient 
period. 



zero. This will lead to large energy gradients in the early timesteps of the simulation, and the 
network needs some time (which is proportional to the wall thickness) to wash away these initial 
conditions. Our simulations start at a conformal time 770 = 1 and evolve in timesteps A77 = O.25770 
until a conformal time equal to half the box size (that is, 77 = 1024). 

The conformal time evolution of the co-moving correlation length of the network £. (specif- 
ically A/V cc , A being the comoving area of the walls) and the wall velocities (specifically 
yv, where y is the Lorentz factor) are directly measured from the simulations, using techniques 
previously described in II 1211 . However here we use a newly parallelized version of the code, 
optimized for the Altix UV1000 architecture of the COSMOS Consortium's supercomputer. 



3. Analytic model 

In order to model a defect network one starts from the microscopic equations of motion (the 
Nambu-Goto equations, in the case of strings) and, through a suitable averaging, arrives at 'ther- 
modynamic' evolution equations. The non-trivial part of this procedure is the inclusion of terms 
to account for defect interactions and energy losses. Such terms must be added in a phenomeno- 
logical way, and for their calibration one must resort to numerical simulations. 

For cosmic strings, this procedure leads to the velocity-dependent one-scale (VOS) model 
1711 . which has been thoroughly tested against simulations. One can follow an analogous 



procedure both for the case of monopoles [19] and for domain walls. This latter case was first 
studied in lfl2ll . and more recently provided a preliminary calibration; here we will provide a 
more quantitative one. 

The evolution equation for the characteristic wall lengthscale L (which is related to the wall 
density p w via L = cr/p w , where <x is the domain wall energy per unit area) and their RMS 
velocity v, are as follows 

= (1 +3v 2 )HL + c w v (4) 

at 

d 4= a -^-3Hv). (5) 



dt \L 
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Here c w and k w are the free parameters: the former quantifies energy losses, while the latter 
quantifies the (curvature-related) forces acting on the walls. To a first approximation, these are 
expected to be constant. Note that in the context of the VOS model the characteristic length scale 
L can further be identified with the physical correlation length £ p / n J . The comoving version of 
this was defined in the previous section, and the two are related via 

iphys = a^ c , (6) 
and we are therefore assuming that = L. Note that if 

fcKT, 1 - 5 , (7) 

then, for an expansion rate A defined as before, 

fphys oc t-* 1 -* . (8) 

Neglecting the effect of the wall energy density on the background (specifically, on the Fried- 
mann equations) — the relevant case for our numerical simulations — one can show that the attrac- 
tor solution to the evolution equations ( 1415b is a linear scaling solution 

L — et , v — const . (9) 

Assuming that the scale factor behaves as a oc t A , the linear scaling constants above take the 
following detailed form: 

6 ~ 3A(l-A) ( ' 

v 2 = ~^T~T — ■ (ID 

3/t k w + c w 

Numerically, we look for the best fit to the power laws 

A 1 

— azp w az — Qcrf, (12) 

' be 

yv ocrf ; (13) 

for a scale-invariant behavior, we should have fi = — 1 and v = 0. The dynamics at the beginning 
of the simulation will be dominated by the initial conditions, while for the walls to be sufficiently 
well defined (which is important for accurately measuring walls areas and velocities) the co- 
moving correlation length should be significantly larger than the wall thickness. Our choice of 
the reliable period for the fits is done by inspection of each set of simulations, using these criteria 
O]. Since we end all the simulations when the horizon becomes half the box size, the periodic 
boundary conditions should have no influence on our results. 

In addition to measuring the scaling exponents /i and v, we will also be interested in the asymp- 
totic values of ^ c /j] and yv, which can be related to the macroscopic parameters of the analytic 
model. These are calculated from the last few timesteps of each simulation, on the assumption 
that by then the network has reached scaling — note that our measured scaling exponents are con- 
sistent with this assumption. These are then used to calibrate the model. The scaling solution is 
parametrized by the expansion rate A (such that a oc t A ) and the phenomenological parameters c w 
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Figure 2: The distribution of the clock time needed to run each of our 90 2048 3 simulations on COSMOS, using 256 
CPUs. The time of each ran has been rounded off to the nearest ten-minute bin. 



and k w . Given these parameters, the predicted values for e and v are given by Eqs. ( [TOIfTTI i. and 
therefore we can trivially obtain the value of v, while e is given by 

€=- = — ; (14) 

from these one finally obtains the numerically measured values of c w and k w . 



4. Simulation results and model calibration 

We have carried out three series of 2048 3 simulations in the radiation era (corresponding to 
A = 1/2), the matter era (A = 2/3) and a fast-expansion era with A = 4/5. Each series consisted 
of 30 different simulations with different (random) initial conditions. These were run on the 
COSMOS supercomputer, using 256 CPU, and each took about 7 hours of clock time. The 
distribution of these times is shown in Fig. [2] The distribution of simulation times is mainly the 
result of the fact that the CPUs involved may not be physically close together on the machine: a 
simulation will start when (among other things) any 256 CPUs are available, and the closer they 
are to each other the faster the run time. (The time spent in identifying the walls can also vary 
from run to run. but we believe that this is a smaller effect.) 

Table [T] and Fig. [3] show the compared results of our series of runs. As can be seen from 
the values of the exponents fi and v, we find no deviation from the theoretically expected scaling 
solution. The behavior of the scaling properties in each of the three epochs is also to be expected. 
Smaller Hubble damping (corresponding to a smaller value of A) allows the walls to have higher 
velocities. Therefore they will have more interactions in a given time, which leads to more energy 
losses, and consequently a smaller density (or equivalently a larger correlation length). 

The oscillations in the wall velocities ar early times in Fig. [3]are also worthy of notice. These 
arise as the network is relaxing away from the initial conditions and approaching the scaling 
solution. The fact that is is clearly visible here is in part due to the larger volume of the simula- 
tions (which provides better statistics), but mostly to the fact that we do not resort to an artificial 
damping period to make this relaxation faster. 
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Figure 3: The evolution of the domain wall density (p w = A IV) and velocity ((yv) 1 ) as a function of conformal time, for 
the radiation era (A = 1/2), the matter era (A = 2/3) and the fast expansion era (A = 4/5). The plotted curves are the 
average of 30 different simulations with different (random) initial conditions. 



A 


Fit range 




V 




yv 


1/2 


126-626 


-1.00 + 0.04 


-0.00007 + 0.00004 


0.57 + 0.01 


0.45 + 0.05 


2/3 


31-626 


-0.99 + 0.02 


-0.00003 + 0.00002 


0.51+0.01 


0.37 + 0.03 


4/5 


63.5-751 


-1.01+0.03 


-0.00000 + 0.00001 


0.42 + 0.01 


0.29 + 0.02 



Table 1: The measured scaling properties of 2048 3 PRS numerical simulations of domain wall networks with different 
expansion rates (parametrized by A). The scaling exponents /j and v were fitted using only the latter part of each sim- 
ulation; the fit range (in conformal time) used for each set is shown in the second column. (See 1 6] for a discussion of 
the choice of this range.) The last two columns show the directly measured asymptotic values of g c /Tj and yv, which 
can be related to the macroscopic parameters of the analytic model. All error bars are one-sigma statistical uncertainties 
obtained by averaging over 30 simulations with different random initial conditions. 




Figure 4: The likelihood contours for the VOS model parameters c w and k w , using the directly measured asymptotic 
values of t; c /r] and yv shown in Table[T] 



6 



The measured values of if c /r] and yv in the three epochs can now be jointly used to calibrate 
the VOS model. Through a standard chi-squared minimization procedure we can easily obtain 
the likelihood contours in the c w - k w plane shown in Fig. [4] and through marginalization we 
then obtain our final calibrated model parameters 

c w = 0.34 ±0.16 (15) 

and 

k w = 0.98 ± 0.07 ; (16) 

these results are remarkably consistent with an early qualitative analysis in [12], which had found 
c w ~ 0.5 and k w ~ 0.9, and with our own estimate in [6], which used sets of 1024 3 runs and found 
c w = 0.5 ±0.2 and k w = 1.1 +0.3. 

The fact that the error bars on the model parameters are somewhat larger than those for ^ c /t] 
and yv are an indication of the fact that taking c w and k w as constants is only an approximation, 
though a reasonable one in the context of a one-scale model. In other words, doing the fit sep- 
arately for the three epochs would lead to slightly different best-fit parameters — cf. Table IV in 
[@]. In particular, we might expect that k w is velocity-dependent, as was found to be the case for 



cosmic strings [17]. 



5. Conclusions 

We have carried out the largest field theory simulations of standard domain wall networks 
undertaken to date, with simulation boxes of size 2048 3 , and confirmed that a scale-invariant 
evolution of the network, with % p h ys x t (or equivalently £ c oc rj) and v = const, is the attractor 
solution. 

By comparing the network's asymptotic scaling properties in several cosmological epochs, we 
obtained an accurate calibration for the two free parameters in the velocity-dependent one-scale 
model for domain walls, one of which describes the network's energy loss rate while the other 
describes the (curvature-related) forces acting on the walls. Our result for these are respectively 
c„ = 0.34+0.16 and k w = 0.98+0.07, which are significantly more precise than earlier estimates. 

This combination of analytical and numerical techniques has also been used for studying cos- 
mic strings il8ll2011 and semilocal strings [21 ], although in both cases the model calibration relies 
on significantly smaller simulations than the ones in the present paper. The larger number of de- 
grees of freedom (corresponding to additional scalar fields, although for cosmic strings one has 
the alternative of Goto-Nambu simulations) makes the study of these models computationally 
challenging since in most cases the factor limiting the size of the boxes that can be simulated is 
memory rather than time, but otherwise our methods are directly applicable there. 

As previously said, domain walls are merely the simplest defect one can simulate. There 
remain certain types of well-motivated string (or hybrid) objects for which their physics, and 
most notably their cosmological evolution, remains relatively unexplored. An example is the 
evolution of cosmic strings carrying currents 11221 . even though such strings have been predicted 
in situations where more than one cosmological phase transition is considered, and they also 
arise naturally in SUSY models that predict strings. String currents also appear in the effective 
4D description of higher dimensional strings, as the string position in the internal dimensions 
is described by worldsheet scalar fields giving rise to currents. Finally, there are superstrings 
models (or toy models thereof) where there are junctions connecting string segments which have 
different masses per unit length and different intercommuting properties. 
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All of these more realistic models have more degrees of freedom than domain walls, and 
simulating them therefore requires, for a given box size, an amount of memory that increases at 
least as the number of degrees of freedom. Additionally, the extra complexity also implies that 
each evolution timestep will take longer in these models. Nevertheless, available codes seem to 
scale reasonably well, at least in the range of conditions probed so far. There is therefore a case 
for carrying out simulations using ca. 1000 cpus and 1Tb of memory. 
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